19  Advanced Spatial Analysis

Today we will be focusing on the practice of fancy geospatial data analysis and visualization.

19.1 Overview

Today will cover three things.

  • Layer Controls
  • Buffers

19.1.1 Libraries

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.3     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

19.2 Data

We’ll use the CalEnviroScreen dataset (SoCalEJ) from other classes and the warehouses dataset for San Bernardino for this example.

Let’s load the old SoCalEJ but only the San Bernardino County census tracts. We’ll call it SBDCoEJ to differentiate it.

URL.path <- 'https://raw.githubusercontent.com/RadicalResearchLLC/EDVcourse/main/CalEJ4/CalEJ.geoJSON'
SBDCoEJ <- st_read(URL.path) |>  
  filter(County == 'San Bernardino') |> 
  st_transform("+proj=longlat +ellps=WGS84 +datum=WGS84")
Reading layer `CalEJ' from data source 
  `https://raw.githubusercontent.com/RadicalResearchLLC/EDVcourse/main/CalEJ4/CalEJ.geoJSON' 
  using driver `GeoJSON'
Simple feature collection with 3747 features and 66 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 97418.38 ymin: -577885.1 xmax: 539719.6 ymax: -236300
Projected CRS: NAD83 / California Albers

Now let’s load the SBD County warehouses.

WH.url <- 'https://raw.githubusercontent.com/RadicalResearchLLC/WarehouseMap/main/WarehouseCITY/geoJSON/comboFinal.geojson'
SBDwarehouses <- st_read(WH.url) |>  
  filter(county == 'San Bernardino') |>  
  filter(category == 'Existing') |> 
  st_transform("+proj=longlat +ellps=WGS84 +datum=WGS84")
Reading layer `comboFinal' from data source 
  `https://raw.githubusercontent.com/RadicalResearchLLC/WarehouseMap/main/WarehouseCITY/geoJSON/comboFinal.geojson' 
  using driver `GeoJSON'
Simple feature collection with 9295 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -118.8037 ymin: 33.43325 xmax: -114.4085 ymax: 35.55527
Geodetic CRS:  WGS 84

19.2.1 Example 1 - Layers Control on Leaflet maps

Let’s make a leaflet map showing the distribution of Hispanic census tracts in San Bernardino County.

First, we need to define the color palette. Feel free to choose a different palette. I’ll choose five quantile bins to get a round 20% number.

palHispanic <- colorQuantile(palette = 'BuPu',
                             domain = SBDCoEJ$Hispanic,
                             n = 5)

Figure 19.1 shows the distribution of census tracts with a legend in a leaflet map.

leaflet()  |>  
  addProviderTiles(providers$CartoDB.Positron) |> 
   setView(lat = 34.1, lng = -117.34, zoom = 10) |> 
  addPolygons(data = SBDCoEJ,
              color = ~palHispanic(Hispanic),
              weight = 1) |> 
  addLegend(data = SBDCoEJ,
            pal = palHispanic,
            title = "% Hispanic",
            values = ~Hispanic)
Figure 19.1: Percentage of Hispanic residents by census tract

Pretty good.

Now let’s make a map with warehouses. We’ll make the warehouses black and show the same area.

Figure 19.2 shows the result.

leaflet()  |>  
  addProviderTiles(providers$CartoDB.Positron) |> 
   setView(lat = 34.1, lng = -117.34, zoom = 10) |> 
  addPolygons(data = SBDwarehouses,
              color = 'black',
              weight = 1,
              fillOpacity = 0.7) |> 
  addLegend(data = SBDwarehouses,
            colors = 'black',
            labels = 'Existing warehouses')
Figure 19.2: Warehouses in SBD County

Now let’s combine the two maps and show both at the same time.

All we need to do here is pull the AddPolygons() and AddLegend() functions from both individual maps into a single code chunk. No new code at this point. Note that each call has to include the data we are pointing to individually.

Figure 19.3 shows the combined overlays.

leaflet()  |>  
  addProviderTiles(providers$CartoDB.Positron) |> 
   setView(lat = 34.1, lng = -117.34, zoom = 10) |> 
  addPolygons(data = SBDwarehouses,
              color = 'black',
              weight = 1,
              fillOpacity = 0.7) |> 
  addLegend(data = SBDwarehouses,
            colors = 'black',
            labels = 'Existing warehouses') |> 
  addPolygons(data = SBDCoEJ,
              color = ~palHispanic(Hispanic),
              weight = 1) |> 
  addLegend(data = SBDCoEJ,
            pal = palHispanic,
            title = "% Hispanic",
            values = ~Hispanic)
Figure 19.3: Warehouses and percent Hispanic

Now let’s add in the layers control to toggle the individual layers on and off.

This requires one new function addLayersControl(). It includes an overlay option called overlayGroups. Lastly we add a new argument in the existing code chunk that defines each layers as a group. Note that including group identifiers for both Polygons and Legends is required to have them both toggle on and off.

Figure 19.4 shows the working toggles as a little stack in the upper right corner above the legends.

leaflet()  |>  
  addProviderTiles(providers$CartoDB.Positron) |> 
  setView(lat = 34.1, lng = -117.34, zoom = 10) |> 
  addLayersControl(overlayGroups = c('Warehouses', '% Hispanic')) |> 
  addPolygons(data = SBDwarehouses,
              color = 'black',
              weight = 1,
              fillOpacity = 0.7,
              group = 'Warehouses') |> 
  addLegend(data = SBDwarehouses,
            colors = 'black',
            labels = 'Existing warehouses',
            group = 'Warehouses') |> 
  addPolygons(data = SBDCoEJ,
              color = ~palHispanic(Hispanic),
              weight = 1,
              group = '% Hispanic') |> 
  addLegend(data = SBDCoEJ,
            pal = palHispanic,
            title = "% Hispanic",
            values = ~Hispanic,
            group = '% Hispanic') 
Figure 19.4: Warehouses and percent Hispanic with overlay toggles

Voila!

19.2.2 Example 2 - Create a buffer around a point or polygon.

The function st_buffer() from the sf package allows one to generate larger shapes surrounding an existing point, line, or polygon.

It only requires the existing dataset and a distance. The only tricky part is figuring out the unit one is buffering in. In some coordinate reference systems, you are in meters. In others, you are in decimal degrees. The buffer distance of 1000 meters is useful, but 1000 latitude degrees goes off map.

Allow me to demonstrate a buffer map with the Gardens dataset. We’ll do a coordinate transformation to NAD83 (a meter based projection) and then back.

gardens <- st_read(dsn = 'CommunityGardens.geojson') |> 
  filter(County == 'San Bernardino') |> 
  #transform to NAD83 meters based coordinates
  st_transform(crs = 4269)
Reading layer `CommunityGardens' from data source 
  `C:\Dev\EA078_Fall2023\CommunityGardens.geojson' using driver `GeoJSON'
Simple feature collection with 186 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -119.2673 ymin: 32.78017 xmax: -115.5097 ymax: 34.52903
Geodetic CRS:  WGS 84
buffer_gardens <- gardens |> 
  st_buffer(dist = 1000) |> 
  #transform back to WGS84 for map
  st_transform(crs = 4326)

leaflet() |> 
  addTiles() |> 
    setView(lat = 34.1, lng = -117.34, zoom = 10) |> 
    addPolygons(data = SBDCoEJ,
              color = ~palHispanic(Hispanic),
              weight = 1,
              group = '% Hispanic') |> 
    addPolygons(data = buffer_gardens,
                color = 'darkgreen',
                fillOpacity = 0.8,
                label = ~Name)
Figure 19.5: San Bernardino County census tract 1000 m buffers

19.2.2.1 Class Exercise

  • Add a layersControl() to this map to toggle on and off the gardens and Hispanic layers.